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Abstract: 

2-dimensional  and  1-dimensional  computer  simulations  of  shock-wave 
phenomena  and  uniform  high-rate  deformation  in  composite,  porous 
and  brittle  materials  have  been  performed.  A  resonance  regime  of 
shock  compression  of  1-D  laminated  model  of  the  composite  is 
discussed.  An  empirical  constitutive  relationship  has  been 
constructed  to  describe  the  dispersive  effects  in  a  composite  in 
continuous  manner.  Computer  simulation  of  dynamic  compaction  shows 
changing  symmetry  of  pores.  The  formation  of  shear  bands  of 
various  configurations  with  friction  and  without  friction  was 
analyzed  in  two-dimensional  geometry.  It  has  been  shown  the  shear 
bands  can  be  responsible  for  the  multi-wave  structure  forming  at 
shock  compression  of  brittle  material.  The  constitutive  model  of 
brittle  material  has  been  proposed  and  tested  in  series  of  1-D 
hydrodynamic  calculations.  The  model  includes  two  parallel 
elements.  One  of  them  works  as  an  usual  elastic-plastic  body, 
another  element  describes  the  friction  and  dilatancy  in  the 
comminuted  component.  The  comminuted  component  fraction  is  a 
scalar  damage  variable  which  grows  from  a  value  of  0  to  1  as  a 
result  of  inelastic  deformation  of  the  intact  component. 
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1.  INTRODUCTION 


Ceramics,  composite  and  porous  materials  are  widely  used  as 
structural  materials  in  the  space  and  rocket  technology.  To 
predict  consequences  of  intense  attacks  like  micrometeorit 
impacts,  it  is  necessary  to  describe  the  behavior  of  similar 
materials  under  shock-wave  loading.  Usual  way  to  analyze  such 
phenomena  is  its  computer  simulation  where  materials  properties 
are  described  by  equations  of  state  and  constitutive 
relationships.  Current  hydrocode  modeling  capabilities  are 
confined  almost  exclusively  to  homogeneous  isotropic  media  and 
deal  with  a  macroscopic  continuum.  It  seems  that  constitutive 
relationships  for  routine  using  in  hydrocodes  have  to  be  continual 
also. 

For  heterogeneous  materials,  which  consist  of  two  or  more  phases, 
several  factors  complicate  constitutive  relation  constructing. 
Response  of  such  kind  of  materials  is  determined  by  behavior  of 
individual  components  as  well  as  their  interactions  and  has  to 
depend  on  the  space  distribution  of  components  in  the  material.  A 
characterization  of  the  dynamic  behavior  of  heterogeneous  material 
requires  more  profound  understanding  of  the  phenomena.  Current 
hydrocode  modeling  capabilities  are  confined  almost  exclusively  to 
homogeneous  isotropic  media  and  deal  with  a  macroscopic  continuum. 
It  seems  that  final  constitutive  relationships  for  routine  using 
in  hydrocodes  have  to  be  continual  also.  On  the  other  hand,  a 
computer  simulation  of  processes  in  model  heterogeneous  materials 
can  give  a  necessary  information  to  construct  these  relationships. 

The  main  objective  of  this  work  is  the  analysis  and  description  of 
the  behavior  of  composite,  porous  and  brittle  materials  under 
impact  loading.  The  goal  is  to  construct  models  of  heterogeneous 
materials  and  constitutive  relationships  for  computer  simulation 
of  impact  phenomena.  The  work  includes  two-dimensional  computer 
simulation  of  deformation  processes  in  models  of  these 
heterogeneous  materials.  Two-dimensional  computer  simulation  gives 
a  detailed  picture  of  the  process  on  the  microscopic  level. 
Results  of  the  simulations  are  then  used  to  construct  constitutive 
relationships  that  are  incorporated  then  into  the  hydrodynamic 
codes . 
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The  Part  2  of  this  report  presents  the  state  of  the  problem.  It 
includes  a  literature  overview  of  experimental  and  theoretical 
results  in  field  of  response  of  composite,  porous  and  brittle 
materials  to  the  impact  loading  and  some  analysis  of 
compressibility  of  mixtures.  Part  3  presents  results  of  2-D 
computer  experiments.  Part  4  contains  constitutive  models  and  data 
of  1— D  calculations.  Mathematical  models  and  codes  are  described 
in  Appendix. 


2.  BACKGROUND  AND  LITERATURE  OVERVIEW 
2.1.  COMPOSITE  MATERIALS  CONTINUUM  MODELS. 


Series  of  simplified  approaches  have  been  developed  to  describe  a 
viscous- elastic  behavior  of  composite  materials  at  low  and 
moderate  strain  rates  [1-4].  These  approaches  are  based  on  various 
methods  of  the  mechanical  parameters  averaging.  For  example, 
effective  modules  theory  [1,2]  uses  assumption  that  average 
stresses  and  average  strains  8^  are  related  by  the  Hooke's 
low  through  the  effective  modules  C  x 


Effective  modules  depend  on  elastic  properties  of  components  as 
well  as  on  the  geometrical  form  and  space  distribution  of 
components.  This  theory  is  based  on  an  approach  that  all 
components  of  stress  tensor  and  strain  tensor  averaged  in  some 
volume  are  equal  to  corresponding  tensor  components  for  composite 
material  taken  as  a  whole.  Effective  toughness  theory  and  its 
modifications  are  used  for  description  of  composites  with  fibers. 
It  is  assumed  the  f ibers-matrix  interaction  occurs  due  to  matched 
displacement  of  fiber  axis  and  matrix.  Kinetic  and  potential 
energies  of  the  composite  deformation  are  the  sum  of  corresponding 
energies  of  the  composite  components,  which  are  taken  with  weights 
assigned  in  according  to  volume  fraction  of  components. 
Anisotropic  model  of  composite  materials  involving  strength 
effects  is  discussed  in  ref.  [5]. 

In  conditions  of  shock-wave  loading  materials  work  mainly  in  the 
plastic  deformation  region.  As  a  result,  many  factors  considered 
in  theories  mentioned  above  have  an  insignificant  influence  on  the 
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dynamics  of  media.  On  the  other  hand,  these  theories  are  based,  as 
a  rule,  on  the  acoustic  approach  and  can  not  be  applied  for 
description  of  intense  dynamic  loading.  An  additional  factor  of 
the  impact  loading  of  composites  is  the  acoustic  interaction 
between  matrix  and  fibers.  Wave  reflections  between  components  of 
the  material  occur  as  a  result  of  their  different  dynamic 
impedances.  Wave  reverberations  produce  a  dispersion  of  a  load 
pulse  in  the  matter.  This  leads  to  changing  attenuation  of  the 
load  and  has  to  be  taken  into  account  at  analysis  of  hypervelocity 
impact  phenomena.  In  order  to  calculate  the  composite  material 
response  to  shock— wave  loading,  it  is  necessary  to  have  in 
computer  code  an  eguation  of  state  of  the  material  and 
constitutive  relations  describing  dispersion  properties  of  the 
material . 

2.2.  SHOCK  COMPRESSIBILITY  OF  COMPOSITE  MATERIALS. 

The  shock  compressibility  of  mechanical  mixtures  of  two  materials 
with  known  Hugoniots  is  calculated  as  a  sum  of  the  components 
compressibilities  [6].  According  to  this  additive  approach,  the 
specific  volume  V  of  mixture  under  pressure  P  is 

V  (P)  «  a  Vi  (P)  +  (1  -  a)  V2  (P)  , 

where  a  is  the  mass  fraction  of  first  component  in  the  mixture, 
V  (P)  and  V (P)  are  specific  volumes  of  components  under  pressure 
P.  The  most  consistent  consideration  of  shock  compressibility  of 
two— component  mixture  has  been  done  by  G.E. Duvall  and  S.M. Taylor 
[7].  As  a  rule,  the  thermal  equilibrium  is  not  reached  in  real 
situations.  Fortunately,  wave  propagation  is  not  very  sensitive  to 
thermal  effects  in  the  moderate  pressure  region.  Comparison  of 
additive  calculations  with  experimental  data  [6]  shows  reasonably 
well  agreement  between  measured  and  calculated  Hugoniots  of 
mixture  when  Hugoniots  of  components  were  used  to  calculate  Vi  (P) 
and  V2 (P) .  Then  the  shock  front  velocity  and  other  kinematic 
parameters  of  shock  waves  in  the  mixture  are  determined  through 
the  Hugoniot.  Some  formal  approximate  equation  of  state  can  be 
constructed  also  using  the  P(V)  Hugoniot  of  the  mixture  [8] . 
Calculations  of  compressibility  of  the  mixture  using  the  mixture 
theory  [9]  or  some  mechanical  models  [10]  are  much  more 
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complicated  but  provide  approximately  the  same  results. 

In  frames  of  additive  approach  the  bulk  sound  velocity  in  mixture 
C  is  formally  calculated  through  known  bulk  sound  velocities  in 
components  using  relationship 


dV 

dVi 

dV 

V2 

(1) 

-  = 

a  — 

+  (i-a)  — 2 

=  -  — 

OP 

dP 

dP 

c2 

Thus,  a  square  of  the  initial  bulk  sound  velocity  in  mixture  is 
calculated  as 


a  /p2  C\  +  (i-a)/p2  C2z 

It  is  interesting  to  mention  that  sound  velocity  value  for  mixture 
can  be  less  then  sound  velocity  values  for  both  components.  Figure 
1  shows  the  sound  velocity  of  a  aluminum/tungsten  mixture 
calculated  as  a  function  of  the  tungsten  mass  fraction.  The 
dependence  has  a  minimum  near  0.9  mass  fraction  of  tungsten. 
Another  interesting  point  is  a  relationship  between  a  shock  front 
velocity  and  a  particle  velocity  in  mixtures.  It  is  known  this 
relationship  is,  as  a  rule,  linear  one  for  metals,  alloys  and 
other  individual  matters.  Figure  2  shows  Hugoniots  of  mechanical 
mixtures  of  tungsten/aluminum,  tungsten/ epoxy  and  aluminum/ epoxy 
calculated  for  volume  fraction  of  heavy  component  equal  to  0.375 
using  linear  Hugoniots  of  all  components.  One  can  see  Hugoniots  of 
mixtures  are  not  exactly  linear  but  curvature  is  not  significant. 

2.3.  POROUS  MATERIALS 

The  initial  stage  of  the  compression  of  a  porous  body  is 
characterized  by  elastic  strains.  The  yield  strength  decreases 
with  increasing  initial  porosity  of  the  medium.  Irreversible 
compaction  occurs  above  the  elastic  limit.  The  load  required  to 
produce  a  given  density  is  found  to  increase  with  decreasing 
initial  density  of  the  body.  This  is  explained  by  strain  hardening 
of  grains  in  the  material  during  the  process  of  compaction. 
Removal  of  the  load  is  not  accompanied  by  a  large  change  in 
porosity.  Complete  compaction  of  porous  material  is  attained  at 
compression  stress,  which  is  somewhat  higher  than  the  Hugoniot 
elastic  limit  for  solid  material. 
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Several  constitutive  models  were  developed  to  describe  a  behavior 
of  the  porous  media  through  the  solid  materials  properties  and 
some  mechanics  of  a  shock  wave  interaction  with  pores.  Each  model 
contains  several  free  parameters  or  material  constants  that  can  be 
determined  only  by  fitting  these  models  to  experimental  Hugoniot 
data.  The  compaction  models  have  frequently  been  formulated  by 
solving  the  stress-strain  relationships  in  a  specific  porous 
geometry.  A  theoretical  prediction  of  the  unique  equilibrium 
compaction  relation  requires  detailed  knowledge  on  material 
property,  pore  collapsing  mechanism,  and  microscopic  stacking 
geometry  of  matrix  material. 

One  from  first  models  of  porous  body  is  so-called  p-a  model  of 
W. Herrmann  [10].  Parameter  of  porosity  a  in  this  model  is  the 
ratio  of  the  pores  specific  volume  V  to  the  specific  volume  of 
solid  matrix  V  .  A  total  volume  of  the  porous  mater  under  stress 

m 

is  changing  due  to  closing  pores  and  due  to  changing  volume  of  the 
matrix.  The  Herrmann's  relationship  for  porosity  is 

a  =  1  +  «Xy  -  1)  [(ps  -  P)/(PS  "  Yp)f  ' 

where  a^  is  a  porosity  value  at  the  yield  stress,  Ps  is  the 
pressure  of  complete  compaction,  Yp  is  yield  stress  for  the 
incipient  porous  material.  A  thermodynamic  equation  of  state  for 
porous  material  p  =  p(V/a,  E)  is  the  same  like  one  for  the  solid 
material . 

R.R.Boade  [11]  has  found  that  the  exponential  p-a  model  fits 
experimental  data  more  accurately.  The  exponential  p-a  model  was 
given  by 

a  =  l  +  (ay  -  i 

where  a  and  p^  are  the  fit  constants. 

M.  Carroll  and  A.  Holt  [12]  used  an  equation  of  state  in  the  form  p 
=  a1f(V/CL,  E)  which  takes  into  account  decreasing  area  in  the 
body  cross-section  of  the  pressure  acting  due  to  the  porosity. 

A  simple  constitutive  model  for  the  shock  Hugoniot  of  porous 
materials  in  the  incomplete  compaction  regime  was  developed  in 
ref.  13.  The  model  is  based  on  the  assumptions  that:  (1)  the 
compacted  powder  geometry  in  the  incomplete  compaction  regime  is 
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approximately  similar  for  the  shock  compaction  and  the  quasi¬ 
isostatic  compaction,  and  (2)  the  net  effect  of  the  strain  rate 
and  the  high  temperature  developed  in  the  shock  compaction  results 
in  a  material  strength  whose  magnitude  is  similar  to  that 
appearing  in  the  case  of  the  quasi—isostatic  compaction.  With 
these  assumptions  the  constitutive  relationship  has  been  obtained 
in  the  following  form: 

7  -  7  (p)[l  +  (  — i 

o 

where  V  (p)  is  a  specific  volume  of  solid  matrix  as  a  function  of 

pressure,  V  and  Vq  are  initial  specific  volume  values  for  porous 

and  solid  materials,  n  «  2.1  to  2.5  is  the  Meyer  work-hardening 

index  of  solid  material  defined  from  hardness  measurements  at 

different  loads,  0  is  yield  strength  of  the  porous  material. 

y 

Models  described  above  do  not  consider  a  stress  relaxation  during 
compacting.  In  the  ref.  14  by  Carroll  and  Holt  some  simple 
kinetics  of  the  spherical  pores  collapsing  was  incorporated  into 
the  model  using  the  Maxwell's  viscosity: 

.  P-P0(a> 

a  =  1(a)  I  p  +  - - - J , 

where  M( CL)  is  a  function  describing  an  elastic  behavior  of  CL, 

p  (a)  is  an  equilibrium  pressure  for  given  a  value.  Ansamble  of 

spheres  with  radius  i)  each  of  them  contained  an  internal  void 

with  radius  CL  was  considered.  The  porosity  a  is  related  with  CL  , 
o 

t>  through  the  relationship 

<W3  =  v(ao_1)- 

According  to  this  model ,  the  relaxation  time  %  is  proportional  to 
the  pore  size  aQ : 

%  -  [par;  (3Y«xo-  i)2/3]1/2, 

where  p  is  the  solid  matrix  density,  Y  is  yield  stress  of  the 
matrix.  Results  of  computer  simulation  of  the  porous  aluminum  with 
this  model  and  work-hardening  incorporated  to  it  (ref.  15)  show 
generally  good  agreement  between  theory  and  experiment,  including 
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both  quasistatic  and  shock  data.  It  has  been  shown  the  viscosity 
is  more  important  to  describe  the  shock  front  rise  time  than  size 
of  pores.  An  adequate  description  of  the  width  of  a  steady 
compression  wave  in  porous  aluminum  was  achieved  by  using 
viscosities  in  the  range  1-10  Pa  s. 

Model  of  Carroll  and  Holt  was  developed  in  ref.  16  with  more 
careful  description  of  viscosity  and  dynamics  of  pores  collapse. 
Two  material  parameters  K  =  (1/D) (Y/po )1/2  and  R  =  T)/ [dQ  (7po ) 1  /2  ] 
allows  to  analyze  the  behavior  of  mater  over  a  wide  range  of 
stresses  and  material  properties. 

Following  development  of  model  has  been  done  in  ref.  17  through 
the  incorporation  of  deviatoric  stresses  into  the  model.  It  has 
been  done  by  such  a  way  that  porosity  and  evolution  of  the  pores 
form  is  determined  not  only  by  pressure  but  by  full  tensor  of 
stresses.  According  to  model,  pores  can  be  closed  even  by  shear 
stresses  only.  The  constitutive  relationship  is  based  on  the 
conception  of  the  yield  surface  in  the  space  of  stresses  which  is 
varied  as  a  function  of  the  porosity  OU 

The  shock  compression  of  porous  media  produces  greater  heating  of 
the  medium  than  the  shock  compression  of  solid  materials.  The 
distribution  of  dissipated  shock  energy  in  the  medium  depends 
mainly  on  how  the  pores  are  collapsed  during  compression.  The 
porous-body  model,  simulated  by  a  set  of  steel  balls,  has  been 
used  in  ref.  18  to  show  that  the  deformation  produced  during  the 
impact  compression  of  a  porous  body  is  localized  mostly  near  the 
surface  of  granules.  The  theoretical  analysis  of  viscoplastic 
heating  of  matter  in  the  neighborhood  of  collapsing  spherical 
pores  was  carried  out  in  refs.  19  to  22.  A  micro-level  numerical 
simulation  shows  the  overheating  matter  near  the  pores  at  shock 
compression  also.  It  was  demonstrated  in  calculations  of  Mader 
[23]  of  shock-wave  propagation  in  liquid  media  containing  closed 
cavities.  The  shock-wave  interaction  with  density  discontinuities 
leads  to  formation  of  regions  of  high  temperatures.  The  size  of 
similar  "hot  spots"  is  close  to  initial  diameter  of  pores.  The 
same  effect  was  observed  in  ref.  24  where  consolidation  of 
composite  material  and  evolution  of  the  particles  form  at  shock 
compression  were  analyzed. 
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2.4.  BRITTLE  MATERIALS 


The  high-strain-rate  behavior  of  ceramics  and  other  brittle 
materials  has  more  individual  peculiarities  than  do  metals. 
Different  modes  of  a  brittle  materials  deforming  under  compressive 
and  shear  stresses  are  illustrated  schematically  in  Figure  3. 
Oxides  and  intermetallic  compounds  have  a  very  high  energy  of 
nucleation  of  dislocations  which  are  responsible  for  the  plastic 
deformation  mechanisms  of  metals  and  other  ductile  crystalline 
materials.  Instead  of  dislocations,  inelastic  deformation  of 
brittle  materials  is  governed  by  the  cracking.  At  low  pressure  the 
shear  cracks  and  cracks  oriented  parallel  to  the  compression 
direction  nucleate  under  stresses  of  1/3  to  2/3  of  the  elastic 
limit  [25].  The  cracking  is  accompanied  with  some  small  volume 
increment,  —  so  called  dilatancy,  of  order  of  1%  or  less. 

A  number  of  physical  mechanisms  [25]  are  responsible  for  the 
dilatational  response  of  brittle  materials  at  triaxial  compression 
tests.  They  include,  for  example,  the  opening  of  secondary  cracks, 
void  growth  in  the  vicinity  of  the  crack  tips  resulting  from  a 
locally  tensile  stress  field,  and  void  opening  at  the 
intersections  of  cracks.  In  the  case  of  granular  materials, 
compressional  dilatancy  may  result  simply  from  a  loss  of  closed 
compactness.  Initially  separate  micro-cracks  are  not  united  and  do 
not  cause  the  body  fracture  on  the  whole.  The  growth  and  collision 
of  cracks  at  following  increasing  shear  stresses  produce  a 
fracture  of  the  material.  If  pressure  is  high  enough,  a  plastic 
deforming  requires  less  energy  than  cracking  and  brittle  materials 
become  ductile. 

A  deformation  of  solids  in  shock  waves  is  one-dimensional  one  and 
shear  stresses  arise  simultaneously  with  increasing  pressure. 
Depending  on  a  relation  between  elastic  modules  and  the  cracking 
threshold  as  a  function  of  pressure,  the  state  of  material  can 
pass  through  the  dilatation  area  or  come  directly  to  the  area  of 
plastic  deformation.  It  seems,  both  these  cases  were  observed  in 
series  of  experiments  with  plane-shock-wave  loading  of  different 
ceramics  performed  during  last  decade. 

Most  advanced  and  informative  measurements  of  shock  compression 
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and  release  in  high-strength  ceramics  were  performed  by  M.E.Kipp 
and  D.E. Grady  [26-28].  The  shock  load/unload  profiles  of  silicon 
carbide,  titanium  diboride,  boron  carbide  and  zirconium  dioxide 
ceramics  subject  to  plate  impact  provide  distinctive  features  of 
the  material  response  that  can  be  used  in  the  formulation  and 
development  of  constitutive  models.  Important  results  were 
obtained  also  by  the  group  of  Z. Rosenberg,  S. Bless,  and  N.S.Brar 
[29-31]  and  D.P.Dandekar  et  al.  [31-33]  for  the  alumina  and 
titanium  diboride,  W.-D.  Winkler  and  A.J.Stilp  [34]  for  alumina, 
titanium  diboride,  silicon  carbide,  and  boron  carbide.  These 
groups  measured  also  the  spall  strength  at  shock  intensities 
below  and  above  the  Hugoniot  elastic  limit  and  deviator  stresses 
in  shock-compressed  ceramics .  The  post-test  examination  of 
recovered  ceramic  samples  was  carried  out  [32,33,35,36]. 

Results  of  experiments  with  ceramics  exhibit  all  possible  load 
ways  shown  in  Fig. 3.  For  example,  stress-strain  diagram  of  the 
silicon  carbide  is  typical  for  elastic-plastic  materials;  shock 
compression  of  the  boron  carbide  ceramics  is  accompanied  with 
cracking  and,  as  a  consequence,  with  decreasing  shear  strength; 
it  seems  the  behavior  of  the  titanium  diboride  is  some 
intermediate  between  former  two  cases  with  crack  nucleation  start 
below  the  elastic  limit  and  obvious  elastic-plastic  following 
behavior. 

The  theories  called  usually  the  damage  continuum  mechanics  are 
widely  used  for  phenomenological  description  of  high-strength 
brittle  material  response  [37-44].  The  state  of  material  element 
(for  simplicity,  in  isothermal  approximation)  is  characterized  by 
e  -  total  strain  tensor,  ep  -  plastic  (or  viscous)  strain  tensor, 
w  -  tensor's  or  scalar  parameter  of  damage.  Being  the  material 
reaction,  the  stress  tensor  is  supposed  to  be  a  function  of  actual 
strains  and  damage,  i.e. 

G  =  G(e,eP,w) 

System  of  constitutes  equations  is  also  including  the  law  of 
plastic  (viscous)  flow 


dep/dt  =  $(G,ep,v) 
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and  an  equation  of  a  damage  evolution 

dw/dt  =  C2(G,ep,w) 

In  case  of  rate- independent  materials  the  functions  $  and  Q 
depends  also  on  the  rate  of  <7  changing.  Relationships  (2),  (3) 
connect  rates  of  the  plastic  strain  and  damage  evolution  with  the 
actual  state  of  material  element.  The  second  group  of  essential 
assumptions  includes  the  suppositions  concerning  to  existence  of 
yield  surface  and  damage  one,  gradientality  (or  nongradientality) 
of  plastic  flow  and  damage  evolution  process  etc.  The 
phenomenological  description  of  damaged  materials  is  connected 
with  definition  of  damage  measure. 

Series  of  attempts  to  construct  reasonable  constitutive  model  for 
description  of  response  of  hard  ceramics  are  known.  F.L.Addessio 
and  J.N. Johnson  [45]  are  developing  a  microphysical ly  based  model 
of  the  brittle  materials  inelastic  deformation  by  cracking,  which 
includes  the  progressive  loss  of  strength  and  degradation  of  shear 
module  as  well  as  the  post-failure  response  of  a  granular  material 
with  friction.  Crack  instability  conditions  and  inelastic  strains 
are  obtained  by  considering  the  response  of  individual  micro¬ 
cracks  to  an  applied  stress  field.  Similar  model  of  D.J. Grove  and 
A.M.Rajendran  [46]  assumes  preexisting  randomly  distributed 
microcracks  in  ceramics,  plastic  flow  above  the  yield  strength, 
degradation  of  elastic  modules  and  strength  due  to  microcracking 
under  both  compression  and  tension.  D. Steinberg  [47]  has  performed 
a  computer  simulation  of  the  shock-wave  processes  in  ceramics  with 
phenomenological  approach  which  was  successfully  used  earlier  for 
metals.  He  described  the  yield  strength  as  a  function  of  the 
pressure,  temperature,  strain  and  strain  rate.  The  description 
includes  also  the  simplest  form  of  the  Cochran-Guinan  Baushinger 
model  [48]  with  degrading  shear  module  at  reverse  deformation.  The 
phenomenological  model  with  the  damage  softening  term  was  used 
also  in  calculations  of  J.F. Davis  et  al.  [49].  The  damage 
evolution  was  assumed  to  be  a  function  of  the  plastic  work  rate  in 
compression  only. 

One  have  to  say  that  good  agreement  between  results  of  computer 
simulation  and  experimental  profiles  takes  place  just  for  rather 
elastic-plastic  materials  like  the  silicon  carbide.  This  does  not 
take  a  place  for  the  materials,  like  boron  carbide,  which  are 
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crushed  under  compression  with  the  shear  component.  It  seems  the 
advanced  model  should  include  both  the  strain  hardening  and 
softening,  Baushinger  effect,  dilatancy  and  a  factor  of  the  strain 
rate. 


3.  RESULTS  OF  2-D  COMPUTER  EXPERIMENTS. 

Computer  simulation  of  the  heterogeneous  materials  response  to 
impact  loading  has  been  carried  out  with  2-D  Lagrangian  code  on 
triangular  grid.  The  mathematical  model  is  described  in  Appendix  A. 

3.1.  2-D  SIMULATION  OF  SHOCK  WAVES  IN  A  COMPOSITE  MATERIAL. 

On  the  first  stage  we  studied  in  the  hydrodynamic  approach  the 
dissipative  properties  of  aluminum  reinforced  by  tungsten  fibers. 
With  this  goal  2-D  computer  simulation  of  steady  shock  waves  in 
the  composite  has  been  performed.  Steady  shock  waves  are  most 
convenient  for  interpretation  and  analysis  because  the  strain 
rate,  the  rate  of  a  longitudinal  stress  increase,  and  deviator 
stresses  are  related  by  simple  dependencies  in  this  case. 

Dimensions  of  triangular  Lagrangian  grid  were  200x40  with  the  5  \MH 
initial  size  of  cells.  The  volume  fraction  of  hexagonally  placed 
tungsten  fibers  was  equal  to  0.375.  The  average  diameter  of 
tungsten  fibers  was  60  | M.  A  cross-section  of  fibers  at  the 
triangular  grid  was  approximated  by  hexagons.  Rectangular  shock 
pulses  of  various  intensities  were  introduced  into  the  material 
through  the  left  boundary  by  linear  increasing  its  velocity  or 
pressure.  Due  to  limited  computer  resources,  it  is  practically 
impossible  to  calculate  with  high  resolution  the  steady  wave 
buildup  from  the  initial  shock  jump.  In  order  to  have  a 
possibility  of  calculations  with  small  cells  of  the  grid  and  to 
keep  an  accuracy  of  calculations,  we  have  performed  series  of 
calculations  with  varied  rise  time  of  the  incipient  load  pulse. 
Results  of  calculations  were  internal  deformation  of  the  material 
and  longitudinal  stress  profiles  in  different  moments  of  time. 

An  example  of  the  digitized  microstructure  after  shock  compression 
is  shown  in  Fig.  4.  Figure  5  shows  an  example  of  integral  pressure 
(longitudinal  stress)  profiles  which  were  constructed  by  averaging 
pressure  values  in  cells  of  greed,  including  viscous  components. 
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along  series  of  cross-sections  perpendicular  to  the  load 
direction.  The  weight  averaging  was  done  taking  into  account  the 
fraction  of  crossed  cell.  Due  to  limited  amount  of  fibers  in  the 
composite  and  to  finite  grid  dimensions  the  calculated  pressure 
profiles  are  not  monotonous.  It  is  necessary  to  mention  that  in 
case  of  pure  aluminum  the  compression  wave  evolution  was 
monotonous  and  the  calculated  minimal  shock-front  thickness  was 
approximately  equal  to  4  or  5  cell  sizes,  that  is  20-25  [MR.  The 
compression  wave  thickness  is  much  more  in  the  case  of  composite 
materials.  This  means  the  numerical  and  approximation  viscosities 
do  not  have  any  essential  contribution  into  the  digitized  rise 
times  for  the  composite  material  case. 

Calculations  with  the  initial  rise  time  100  ns  provide  the  steady 
compression  wave  buildup  at  1  mm  of  the  propagation  distance.  A 
steadiness  of  the  compression  wave  was  verified  through  its 
average  propagation  velocity.  The  propagation  velocity  of  steady 
wave  equals  to  the  shock  front  velocity  calculated  for  the 
aluminum-tungsten  mixture.  The  thickness  of  steady  compression 
waves  in  composite  amounts  to  several  (3  to  5)  periods  of  fibers 
layers  and  decreases  gradually  with  increasing  the  shock 
intensity. 

Figure  6  shows  stress  profiles  which  were  calculated  taking  into 
account  elastic-plastic  properties  of  both  matrix  and  fibers 
materials.  The  yield  strength  Y0  was  0.5  GPa  for  aluminum  and  3 
GPa  for  tungsten.  Calculations  show  much  more  steep  middle  part  of 
the  compression  wave  in  the  case  of  elastic-plastic  behavior  in 
the  low  pressure  range.  Obviously  there  are  two  reasons  of  the 
fact:  (i)  an  increased  average  sound  velocity  accelerates  the  wave 
reverberations;  (ii)  a  resistance  to  shear  strains  in  solids 
promotes  to  faster  buildup  of  a  final  configuration  of  the  shock- 
compressed  composite.  Only  very  weak  precursor  is  observed  in  the 
case  of  2-dimensional  elastic-plastic  composite,  but  the  final 
stress  exceeds  the  shock  pressure  in  corresponding  liquid-like 
material.  It  seems  the  difference  in  average  final  axial  stresses 
is  controlled  by  the  yield  strength  of  the  matrix. 

3.2.  POROUS  MATERIALS. 

Some  preliminary  results  for  porous  materials  were  obtained  in  two 


13 


modes  of  dynamic  compaction.  2-D  computer  simulation  was  performed 
in  the  same  geometry,  as  shown  in  Figure  4,  with  a  gas  placed 
instead  of  fibers  in  composite  material.  Solid  component  of  porous 
was  described  as  an  ideal  elastic- plastic  body  without 
strain-hardening . 

Figure  7  shows  a  deformation  of  the  grid  in  the  initial  stage  (71 
®^®r  the  process  start)  of  shock  wave  propagation  in  porous 
aluminum.  The  boundary  velocity  was  1  km/s.  Contrary  to  usual 
models  of  shock  compaction,  one  can  see  a  non  spherical  closing  of 
pores  at  one- dimensional  shock  compression. 

A  computer  simulation  permits  one  to  study  a  uniform  compression 
the  material.  For  this  mode  the  uniform  velocity  distribution 
along  the  sample  axis  was  used  as  an  initial  condition  and  the 
constant  values  of  velocity  (uO  for  the  left  boundary  and  0  for 
the  right  boundary)  were  used  as  the  boundary  conditions.  Figure  8 
shows  the  sample  deformation  after  80  ns  at  «  =3  km/s.  In  this 

case  pores  are  closed  more  symmetrically  in  comparison  with  shock- 
wave  compaction,  but  nevertheless  the  spherical  symmetry  is  not 
kept.  Figure  9  shows  the  average  axial  stress  in  porous  tungsten 
as  a  function  of  strain  at  uniform  compression.  After  small 
initial  elastic  part  of  the  diagram  a  remarkable  growth  of 
stresses  is  beginning  just  near  the  end  of  the  compaction  process. 
Obviously  this  is  because  the  strain— hardening  and  viscosity  were 
not  involved  in  calculations. 

The  elastic  limit  of  one-dimensional  compression  is  1.9  GPa. 
Meanwhile  the  yield  strength  of  tungsten  used  in  calculations  was 
3  GPa  what  corresponds  to  5  GPa  of  the  Hugoniot  elastic  limit  of 
solid  tungsten.  According  to  Steinberg's  model  [56]  the  elastic 
limit  of  porous  tungsten  was  expected  to  be  1.70  GPa.  Results  of 
calculations  some  exceed  preestimated  value,  but  the  discrepancy 
is  not  very  large.  It  is  possible  to  say  that  simple  Steinberg's 
model  can  be  quite  used  as  a  constitutive  model  of  porous  elastic- 
plastic  body.  To  take  into  consideration  the  strain  hardening 
effect,  it  is  necessary  to  use  modified  models,  like  model  of  Ki- 
Hwan  Oh  and  Per-Anders  Persson  [ 13 ] . 


3.3.  2-D  SIMULATION  OF  SHOCK  COMPRESSION  OF  A  BRITTLE  MATERIAL 

The  dynamic  compression  of  brittle  material  was  simulated  in  our 
calculations  artificially  as  a  shear  bands  formation  in  the  ideal 
elastic-plastic  body.  Shear  bands  are  formed  in  preliminary 
determined  places  of  the  sample  when  the  stressed  state  in  sample 
reached  some  threshold  conditions  (usually  half  of  the  yield 
strength  of  the  main  material) .  After  that  the  resistance  to  a 
shear  deformation  inside  the  bands  was  set  to  be  zero  or  was 
proportional  to  pressure  like  a  friction.  The  friction  coefficient 
was  considered  as  a  free  parameter  and  was  varied.  At  high 
pressure  the  friction  forces  can  exceed  the  yield  strength.  In  our 
calculations  we  limited  the  friction  resistance  to  shear  by  the 
yield  strength  of  the  ceramic.  An  equation  of  state  and  elastic- 
plastic  parameters  of  the  SiC  ceramic  were  used  in  all  2-D 
calculations . 

Two  series  of  calculations  with  various  geometry  of  shear  bands 
have  been  performed.  An  example  of  deformation  of  a  brittle  sample 
with  intersecting  shear  bands,  is  shown  in  Figure  10.  Figure  11 
shows  the  stress-strain  diagrams  for  samples  with  different 
frictions  inside  shear  bands.  In  general,  these  diagrams  are  quite 
similar  to  observed  experimentally  [26,27]  for  ceramics.  In  the 
case  of  zero  friction,  hydrostatic  states  are  created  inside  the 
shear  bands.  As  a  result,  after  the  shear  bands  formation  the 
structure  of  the  solid  ceramic  is  changed  to  structure  of  solid 
grains  surrounded  by  liquid  layers.  As  a  result,  the  stressed 
state  of  the  material  on  the  whole  becomes  the  hydrostatic  one 
also.  We  did  not  see  any  interactions  between  grains  and 
dilatancy  effects  in  our  model.  In  the  case  of  non-zero  friction 
the  stressed  state  after  the  shear  bands  formation  is  controlled 
by  the  friction  coefficient. 

Besides  the  geometry  shown  in  Fig.  10,  simulations  with  non¬ 
crossing  bands  were  performed  at  the  same  relative  area  of 
inclined  shear  bands.  It  was  expected  that  in  the  last  case  the 
shear  bands  formation  should  not  be  accompanied  with  hydrostatic 
states  of  the  matter.  Figure  12  shows  the  calculated  stress- 
distance  profiles  for  the  ceramic  in  frame  of  the  pure  elastic- 
plastic  model  (without  shear  bands) ,  and  in  frame  of  the  model  of 
elastic-plastic  matter  with  shear  bands.  Of  course,  the  shear 
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bands  formation  lowers  the  Hugoniot  elastic  limit.  In  the  case  of 
non-zero  friction,  the  simulation  gives  three-wave  structure  when 
the  friction  forces  inside  the  shear  bands  exceed  the  yield 
strength.  These  three-wave  structures  are  some  similar  to 
experimental  velocity  profiles  measured  for  the  titanium  diboride 
ceramic.  Qualitatively  stress  profiles  are  not  changed  with 
varying  the  bands  geometry. 

4.  CONSTITUTIVE  MODELS  AND  1-D  COMPUTER  SIMULATION. 

4.1.  SHOCK-WAVE  PROCESSES  IN  LAYERED  COMPOSITES. 

Behavior  of  sophisticated  systems  is  usually  investigated  on 
simple  models.  Such  a  simple  model  of  composite  materials  can  be 
presented  as  a  plate  composed  of  alternating  flat  layers  of  two 
different  materials  perpendicular  to  the  direction  of  shock 
propagation.  Shock-wave  phenomena  in  laminated  composites  were 
discussed  earlier  in  refs.  [50-52].  In  this  case  the  first  shock 
wave  decays  very  rapidly  due  to  numerous  waves  reflections  at 
interfaces.  It  has  been  shown  the  wave  front  induced  a  large 
amount  of  ringings  as  it  passed  through  the  layers  of  the 
composite.  Obviously,  the  final  shock-compressed  state  has  to 
correspond  to  the  Hugoniot  of  the  mixture. 

Figure  13  shows  results  of  computer  simulation  of  shock 
compression  of  laminated  plate  consisting  of  the  copper  and 
polyethylene  layers.  The  pressure  and  particle  velocity 
distributions  at  several  moments  of  time  show  the  resonance 
behavior  of  such  periodical  one-dimensional  composite.  The 
pressure  (and  compression)  oscillations  are  concentrated  mainly 
inside  the  soft  polyethylene  layers.  Due  to  that  one  can  expect 
the  main  dissipation  has  to  be  localized  in  the  soft  material 
also.  The  velocity  oscillations  are  concentrated  mainly  at  more 
heavy  copper  layers.  Figure  14  shows  pressure  histories  in  the 
middle  sections  of  copper-  polyethylene  and  aluminum-polyethylene 
targets.  One  can  see  the  pressure  oscillations  have  an 
approximately  constant  period  and  are  near  harmonic  in  form.  The 
period  of  oscillations  approximately  corresponds  to  the  sound 
propagation  time  through  two  alternating  layers  of  the  composite. 
Wave  reflections  between  layers  smear  the  wave  front,  but  its 
average  propagation  velocity  is  practically  equal  to  the  shock- 
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front  velocity  in  the  mixture  calculated  for  a  given  boundary 
velocity.  An  average  pressure  also  corresponds  to  the  Hugoniot  of 
the  mixture.  The  rise  time  of  first  compression  wave  approximately 
eguals  the  period  of  oscillations.  The  pressure  oscillations  can 
pass  to  a  homogeneous  barrier  placed  behind  of  the  layered  target 
if  the  dynamic  impedance  of  the  barrier  is  high  enough. 

A  reality  of  oscillations  has  been  examined  experimentally.  With 
this  goal  the  pressure  profile  at  the  interface  between  the 
target,  consisting  from  10  copper  foils  0.2  mm  thick  and  10 
polyethylene  films  0.2  mm  thick,  and  copper  plate  was  measured 
with  manganin  pressure  gauge.  The  initial  load  pulse  was 
triangular  in  form  with  a  total  duration  about  20.  Result  of 
measurements  presented  in  Figure  15  confirms  an  appearance  of 
oscillation.  Relatively  small  amplitude  of  measured  oscillations 
and  relatively  fast  decay  of  them  are  explained  by  the  viscosity 
of  soft  polyethylene  layers  where  a  main  deformation  occurs. 

Figure  16  shows  an  influence  of  stress  relaxation  in  the 
polyethylene  layers  on  the  computed  wave  profiles.  One  can  see 
oscillations  decay  much  faster  if  the  relaxation  time  is  some  less 
than  half-period  of  oscillations.  In  the  case  of  large  relaxation 
time  behavior  of  soft  layers  becomes  closer  to  elastic  one,  what 
can  be  seen  as  increasing  the  average  propagation  velocity  and 
decreasing  the  period  of  oscillations.  The  stress  relaxation  in 
heavy  rigid  layers  does  not  produce  any  essential  effect. 

Elastic-plastic  response  of  both  components  increases  the  wave 
dispersion  in  a  composite  and  damps  resonance  oscillations.  This 
can  be  seen  in  Fig.  17  that  shows  calculated  pressure  profiles  at 
the  interface  between  a  tungsten  barrier  and  the  Al/W  layered 
composite  with  the  tungsten  volume  fraction  of  0.375.  The  weak 
precursor  increases  the  total  rise  time  of  the  first  compression 
wave  and  oscillations  decay  much  faster  in  the  elastic-plastic 
composite. 

4.2.  1-D  COMPUTER  SIMULATION  OF  THE  WAVE  DISPERSION  IN  COMPOSITE 
MATERIAL. 

A  dispersion  of  shock  front  takes  place  in  composite  materials  due 
to  waves7  reflections  between  components  with  different  dynamic 
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impedances.  Comparison  of  the  wave  propagation  characteristics  of 
composite  and  viscous  (i.e.,  stress  relaxing)  materials  reveals 
that  there  are  some  basic  similarities  in  their  wave  profiles. 
Composite  and  viscous  materials  appear  to  sustain  the  steady  wave. 
The  smooth  wave  transitions  in  viscous  materials,  which  are 
attributable  to  rate  effects,  appear  to  be  similar  to  smooth  wave 
transitions  in  the  composite  which  are  attributable  to  dispersive 
effects.  In  view  of  this  similarity  the  gross  behavior  of  a 
composite  can  be  described  by  an  appropriate  viscous  material 
constitutive  model,  even  though  the  individual  constituent 
materials  of  the  composite  may  be  strain  rate  independent.  Such  a 
model  was  discussed  earlier  by  L.M. Barker  [51]. 

Following  to  Maxwellian  model,  L.M. Barker  has  introduced  a 
metastable  (instantaneous)  and  equilibrium  stress-strain  paths  of 
loading,  and  stress  relaxation  from  an  instant  state  to  the 
equilibrium  one.  It  is  necessary  to  say  that  determination  of  the 
metastable  loading  path  for  composite  is  a  problem.  On  the  other 
hand,  computations  and  available  measurements  of  shock-wave 
processes  in  composite  materials  do  not  reveal  any  signs  of  the 
ultimate  metastable  states. 

In  this  work  we  tried  to  describe  the  process  with  an  empirical 
constitutive  relationship  without  determination  of  metastable 
paths.  Instead  of  metastable  states  we  assume  that  fast  bulk 
compression  produces  superfluous  pressure  in  matter  which  depends 
on  compression  rate.  Then  an  establishing  the  mechanical 
equilibrium  occurs  through  the  superfluous  pressure  relaxation  to 
zero.  Total  pressure  p  is  presented  as  a  sum  of  equilibrium 
component  pe  determined  by  the  equation  of  state,  and 
nonequilibrium  component  p  which  depends  on  the  strain  rate  and 
load  history.  Nonequilibrium  pressure  components  at  compression 
were  calculated  using  the  empirical  kinetic  relationship 


dP 
_ ] 

dt 


t  dV  )n  P 

_0  _  __  _ n 

V  dt  t 

o 


(C/Co ) , 


(3) 


where  t  =  Is,  K  ,  TL,  and  relaxation  time  T  are  the  material 
o  1 

constants.  An  acceleration  of  the  reverberation  process  due  to 
increasing  the  sound  velocity  with  pressure  is  taken  into  account. 
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The  same  relationship  with  a  reverse  sign  at  (ffl  is  used  for 
rarefaction. 

Empirical  kinetic  relationship  (3)  contains  three  parameters:  Kx  , 
n,  and  T,  ,  which,  in  principle,  can  be  estimated  from  series 
profiles  of  steady  shock  waves  with  different  peak  pressures.  When 
a  steady  wave  propagates  through  a  material,  all  parts  of  the  wave 
propagate  with  the  same  constant  velocity,  and  the  material, 
therefore,  is  loaded  along  a  Rayleigh  line  to  the  final  state  on 
the  equilibrium  behavior  curve.  Intermediate  states  in  the  steady 
compression  wave  contain  some  superfluous  pressures  that  equal  the 
difference  between  stresses  at  the  Rayleigh  line  and  equilibrium 
pressures  corresponding  to  the  material  equation  of  state: 


P„  =  P0  V2  <7,  -  Vi  -  P(V) ;  p o  -  -p2  V  (U2  -  a2), 

dp  1/2 

where  U  is  the  shock  front  velocity,  a  =  (-  V2q  )  is  the 

Lagrangian  sound  velocity.  In  some  intermediate  point  of  the  shock 
wave  the  superfluous  stress  passes  through  the  maximum.  In  this 

point  p  =  0  and,  according  to  (3) , 

L  n 


K  %  =  Fmx  (C/C  ) 

1  n  0 


t  dV  1 
7  dt 

o  ' 


Thus  the  strain  rate  dependence  on  the  p“ax  value  give  an 
estimation  of  the  product  and  the  power  index  71.  A  slope  of 
the  shock-wave  just  in  the  front  point,  where  a  superfluous  stress 
is  steel  zero  and  sound  velocity  is  CQ  ,  gives  us 


to  ^  1 

v  dt 

0 


Last  relationship  can  be  used  to  find  the  K value.  It  shows  also 
the  power  index  71  has  to  exceed  unit  because  in  the  case  of  71=1 
initial  slope  of  the  steady  compression  wave  does  not  depend  on 
its  intensity,  and  in  the  case  of  71  <  1  the  slope  decreases  with 
increasing  shock  intensity. 


Results  of  1-D  simulation  with  Ki  =  1.5'108  Pa/s,  71  =  1.5,  and  T  - 
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2  * 10_S  s  are  compared  with  2-D  data  in  Figures  18.  Equation  of  state 
based  on  the  Hugoniot  of  mixture  was  used  in  1-D  calculation.  One 
can  see  a  sufficient  agreement  between  results  of  two  kinds  of 
calculations. 

4.3.  CONSTITUTIVE  MODEL  FOR  DYNAMIC  RESPONSE  OF  BRITTLE  MATERIALS. 

To  construct  rather  simple  constitutive  model  of  the  brittle 
material  we  used  as  a  basis  the  structural  Marzing  model  [53]. 
This  model  reflects  a  micro-nonuniformity  of  real  materials  and  is 
successfully  used  for  the  Baushinger  effect  description  at  the 
repeated-alternating  loading. 

In  frame  of  the  Marzing  model  model  each  elementary  volume  of  body 
is  represented  as  consisting  from  N  parallel  elastic-viscous- 
plastic  subelements.  Subelements  have  equal  elastic  modulus,  but 
different  yield  strengths  and  viscosities.  The  stress  in  each 
subelement  is  determined  as 

N 

°  -  1  V 
1 

where  0*  are  stresses  inside  the  subelements,  gk  are  weight 
factors  which  can  be  presented  as  relative  areas  of  subelements 
cross-sections.  Inelastic  strain  component  appears  at  moment  when 
an  yield  strength  of  the  weakest  subelement  is  reached.  The 
ultimate  stress  in  this  model  is 

w 

5  -  1  Yk  ek  • 

l 

Effective  elastic  modulus  after  yield  stress  reached  in  n 
subelements  is 

n 

B„.,  -  \  f1  -  l«k). 

l 

where  E  is  the  elastic  modulus  of  subelements.  The  model  can  be 
o 

generalized  on  the  arbitrary  stresses  state  [53]. 

For  high-strain-rate  conditions  it  is  necessary  to  take  into 
account  the  viscosity  of  matter.  Measurements  of  the  shock  front 
structure  [54]  show  the  shear  plastic  strain  rate  in  the  range  of 
104  to  107  s’1  is  related  to  shear  stress  as 
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7  =  a  (t  -  y/2)2, 

p 

where  is  constant.  In  frame  of  the  multi-elements  model  the 

plastic  strain  rate  of  subelements  can  be  calculated  as  [55] 


7 k  =  A  [  -  (T*-  Y.  /2)]2, 

p  y  K 

k 

Both  strain  hardening  and  softening  of  subelements  describe  a 
complex  response  of  materials  to  the  shock-wave  loading. 

Most  of  practical  applications  of  the  computer  simulations  of 
the  brittle  materials  dynamic  response  require  are  associated  with 
three-dimensional  analysis.  From  this  point  of  view  it  is 
necessary  to  avoid  superfluous  details  which  can  be  prohibitively 
expensive  in  realization  of  the  computer  model.  In  relation  to 
multi— elements  model  this  means  that  we  have  to  reduce  to  minimum 
the  number  of  elements.  In  fact,  two  elements  are  quite  enough  to 
describe  with  sufficient  accuracy  the  main  peculiarities  of  the 
material  response. 

Our  model  of  brittle  material  includes  two  parallel  elements.  One 
of  them  works  as  usual  elastic-plastic  body,  possibly  -  with 
strain  hardening.  Another  element  describes  the  resistance  to 
deformation  of  the  comminuted  component.  For  this  second  element 
the  resistance  to  shear  is  a  friction.  Initially  we  deal  with 
solid  undamaged  ceramic,  so  it  is  naturally  to  assume  the  initial 
fraction  of  the  comminuted  component  to  be  zero.  With  beginning  of 
an  inelastic  deforming  the  brittle  material  is  cracking  and  the 
fraction  of  comminuted  component  is  growing  with  increasing 
inelastic  strain.  After  some  strain  the  initial  intact  material 
becomes  completely  failed,  that  means  the  fraction  of  comminuted 
component  reaches  1  and  the  resistance  to  shear  is  determined  by 
friction  in  a  granular  material  only.  Thus  in  frame  of  this  model 
the  comminuted  component  fraction  is  a  scalar  damage  variable 
which  ranges  from  a  value  of  0  to  1.  In  the  model  the  fraction  of 
comminuted  components  grows  with  increasing  inelastic  deformation 
of  intact  component: 

Pf 

r  *P 

B2  “  ~y  f 

2  l  +  671 

r  'p 
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where  71  is  the  inelastic  strain  of  the  intact  component, 
p 

A  dilatational  effect  has  been  formally  incorporated  into  the 
model  to  take  into  consideration  an  increase  in  volumetric  strain 
resulting  from  cracks  growth.  In  the  present  model  the  dilatancy 
is  associated  with  inelastic  deformation  of  comminuted  component 
and  is  expressed  in  terms  of  increments  of  relative  volume  of 
voids.  Obviously  shear  strains  under  compression  can  produce  just 
limited  volume  of  voids.  On  the  other  hand,  positive  pressure 
above  some  threshold  has  to  close  voids.  Basing  on  these  simple 
ideas  we  have  tried  to  describe  the  dilatancy  by  the  following 
empirical  relationship 


c W 


dT2  1  +  ft  V 

1  p  2  v 


(  P  ~  P 

'  p 


where  V  =  V  /  V  is  the  relative  volume  of  voids,  72  and  g  are 
the  inelastic  deformation  and  weight  factor  of  the  comminuted 
component,  ft  ,  ft  ,  p*  ,  PQ  are  the  material  constant  parameters. 
Voids  are  not  closed  when  pressure  does  not  exceed  the  threshold 
value  p* . 


The  model  does  not  include  any  degradation  of  shear  module  with 
inelastic  compression.  The  shear  module  dependence  on  the  pressure 
is  calculated  in  the  approach  of  constant  Poisons  ratio  which 
means  that  ratio  of  the  shear  module  to  the  bulk  module  is  not 
changed  with  changing  pressure  and  temperature.  Yield  strength  of 
the  intact  component  is  changed  proportionally  to  the  shear 
module.  The  tensile  strength  of  the  matter  was  supposed  to  be 
negligible. 


A  computer  simulation  of  shock-wave  experiments  [26,27]  with 
boron  carbide  ceramic  has  been  performed.  The  boron  carbide 
ceramic  was  analyzed  because  previous  calculations  [45-49]  for 
this  ceramic  were  lest  successful.  Results  of  calculations  with 
our  constitutive  model  are  compared  with  experimental  [26,27] 
velocity  profiles  in  Figure  19.  It  seems  the  agreement  between 
results  of  calculations  and  experimental  data  is  quite  reasonable. 
Figure  20  shows  states  of  material  during  shock  compression  and 
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unloading.  One  can  see  the  strong  Baushinger  effect  is  described 
well  by  the  model.  The  dilatancy  effect  is  most  essential  at  last 
stages  of  unloading. 

CONCLUSION 

Results  presented  here  can  be  considered  as  encouraging  in  sense 
that  way  of  investigation  of  the  heterogeneous  materials  response 
to  dynamic  loading  using  the  2-D  computer  simulation  of  the 
phenomena  can  be  really  productive.  Some  important  details  of 
the  dynamic  compression  process  in  composite,  porous  and  brittle 
materials  have  been  established  due  to  2-D  simulation.  Reasonable 
constitutive  relationships  were  then  designed  to  describe  behavior 
of  composites  and  brittle  materials  under  impact  loading.  The 
following  work  can  be  directed  to  more  advanced  analysis  of 
the  deformation  process,  to  development  of  constitutive  models  and 
to  finding  the  models  parameters  values  for  some  determined 
materials. 
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APPENDIX  A. 


MATHEMATICAL  MODEL  FOR  2-D  COMPUTER  SIMULATION 


One  can  hope  that  computer  simulation  of  the  dynamic 
deformation  processes  in  composite  material  can  make  clear  details 
of  these  processes  and  find  correlation  between  mechanical 
properties  of  material  and  its  components.  Some  preliminary 
results  were  obtained  in  our  previous  research. 

Calculations  are  based  on  fundamental  conservation  lows  [1] 
and  are  conducted  in  Lagrangian  material  coordinates.  Conservation 
of  momentum  is  described  by  equations 
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XX  ,  XX 
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dx  +  ~dx~ 


Equation  of  conservation  of  energy  is 


5E 


=  V 
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P  5Y 
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+  V  2  e  +  2  e  + 


XX  XX  zz  zz 


5 

5x 


5T 

dx 


dz 


X 


5T 

dz 


(3) 


The  first  term  in  the  right  part  of  equation  (3)  is  energy  release 
due  to  external  source?  second  term  is  the  work  of  bulk 
strain;  third  term  is  the  work  of  shear  strain;  fourth  term 

describes  the  internal  energy  increment  due  to  thermocondactivity ; 

5U  .  5U  .  5U  5U 

8„=  dx*'  m  '  and  ex,=  dz*  +  dl  are  strain  rates- 

The  mass  conservation  is  carried  out  automatically  due  to 
Lagrangian  variables  used. 


Elastic-plastic  behavior  of  the  material  is  described  by  the 
Hooke's  low  for  strain  rates  [2] 
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A  rotation 

angle  A(p  of  element 

of  media 

in  the  x-z  plane  during 

small  time 

step 

At  is 

au 

au 

Acp  =  At  rot  U  =  -  d^z. 

Corresponding  corrections  in  equations  (4) -(6)  are  calculated  as 
[3] 


6S_  =  (  Szz-  Sxx)sinzA(p  +  2SvvsinAcp  cosA(p 


6S  =  -  6S 


as  =  (  S  -  S  )sinA(p  cosA(p  +  2S  sin2 Acp 

-v  r7  ZZ  XX  Y  V  1 


Yielding  is  tested  using  Mises  criterion.  According  to  this 
criterion,  when 


f  =  2  (  S2  +  S2  +  S2  +  S  S  )  >  §  Y2, 

XX  xz  ZZ  xxzz  o  o 


(T) 


where  Y  is  yield  strength,  yielding  has  occurred  and  the  deviator 
stresses  are  reduced  in  the  usual  way  [81 


S  : 

X  X 
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z  z 
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(8) 

(9) 

(10) 


A  scalar  artificial  viscosity  is  used  to  damp  nonphysical 
oscillations.  Combination  of  the  linear  qx  and  square  qz 
viscosities  is  used  in  our  code: 


<li  =  <li  +  ^  ’ 

where 

qt  =  -  C^AtAZ  V/V  , 


(11  ) 
(12) 
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(13) 


q2  =  CzAZ2  V2/V  , 

hl=y  4  A2/t)  is  a  characteristic  size  of  the  Lagrangian  mesh  [  1  ] ,  A 
is  the  mesh  area,  t>  is  the  sum  of  the  squares  of  the  mesh  sides, 

is  the  dilatation  rate,  Cj-0,5  and  C2^2,0  are 

constant  coefficients.  Artificial  viscosity  is  set  equal  to  zero 
when  dilatation  rate  is  positive. 


V  = 


i  av 

T~5X 


Natural  viscosity  has  been  included  into  the  model  in  form 
[4].  The  bulk  viscosity  is  described  as 


^  +  |r, 


(14) 


q  and  q  values  are  added  to  pressure  in  eqs.  (1)— (3) .  Deviator 
components  of  viscous  stress  are  [4] 
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xX  X 

qzz 
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(15) 

(16) 
(17) 


In  order  to  damp  oscillations  of  stress  components,  deviator 
pseudoviscosities  have  been  introduced  into  the  code  also  as  a 
combination  of  the  linear  [3]  and  square  [1]  regarding  strain 
rates  terms. 


Components  of  the  deviator  of  the  linear  pseudoviscosity 
tensor  are 

q'  =  C,  a  A Z  e  /V,  q’  =  C  a  AZ  e  /V, 

^•xx  3  xx  ^-z  z  3  zz 

q’  =  C  a  AZ  e  /V  , 

^x  z  3  x  z 

and  components  of  the  deviator  of  square  pseudoviscosity  are 

q"  =  C  AZ2  e  H/V  , 

J-Y  V  4  Y  Y 


q"  =  C  AZ2  e*  H/V  , 

^■7  7  4  7  2 


q"  =  C  AZ2  e  ’H/V  , 

^■x  z  4  xz 
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where  H  =  2  /T  is  the  shear  strain  rate  intensity, 

J  =  (6  -  6  )2/6  +  3/2  e2  is  the  second  invariant  of  the 

2  xx  zz  xz 

strain  rates  deviator,  C^O.5,  and  C2^2,0  are  constants.  A  total 
pseudoviscosity  is  sum  of  linear  and  square  components: 

rj 

q  =  q*  +  q"  , 

^X  X  ■Lxx  Axx 

rv 

q  =  qT  +  q"  , 

■L z  z  xz  z  xz  z 

ro 

q  =  q'  +  q"  , 

^x  z  xx  z  ■‘•xz 

Components  of  deviator  of  visco-elastic  stresses  in  eqs.  (l)-(3) 
are  determined  as 


2  =  S  +  q  +  q 

XX  XX  AX  X  Axx 

ro 

2  =  S  +  q  +  q 

z  z  z  z  *z z  xz  z 

ro 

2  =  S  +  q  +  q 

X  z  X  z  ■Lx  z  Xx  z 

An  equation  of  state  of  the  matter  was  used  in  form  of  ref. 
[4]  with  correction  in  field  of  small  densities.  The  pressure  at 


specific  volume 
matter: 


V 


is  calculated  using  Hugoniot  of 


P  =  P„  +  1  <E  "  V  • 


(18) 


where  Ph=  C2(Vo-  V)/CV0-  Sa(VQ-  V)]2  and  Eh=  Ph(Vq-  V)/2  are 
pressure  and  energy  on  the  Hugoniot. 

For  moderate  pressures  in  a  range  of  approximately  0.5  to  50 
GPa,  a  difference  between  the  Hugoniot  and  isentrope  of  the  matter 
is  not  significant  [5],  A  temperature  in  this  region  V  <  VQ  is 

T  =  Th  +  (E  -  Eh)/Cv  , 

where  the  temperature  on  the  Hugoniot  T  is  calculated  as  [8] 

T(v  -  v)  (V  -  V)P„ 

T„  -  T.e  +  °2l!  *  - 
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V 


P  e  TV/V°  [2  (Vo-  V)]  dV 

V 

o 

A  pressure  in  field  of  specific  volumes  more  than  initial 
one  (  V  >  VQ)  is  calculated  using  the  Gruineisen  equation  of 
state 

P  =  ^  (E  -  Es)  ,  (19) 

where 


E  =  (B  K  )/j  (V  -  V) 

s  corr  1  o 

-  ((V/V  )2-  1  ) 

K  =  e 


7  =  2/3+  (7o  -  2/3  )K 


(20) 


The  correcting  factor  K  decreases  the  bulk  modulus  value 

corr 

B  from  B  at  V  =  V  (what  corresponds  to  usual  initial  conditions) 

O  O 

to  zero.  By  this  way  the  energy  of  tension  of  the  condensed  matter 
E  is  reduced  to  zero  at  V  -»  00  affording  a  smooth  transfer  to  the 

S 

gas  state  region.  The  relationship  (20)  for  Kc o r r  is  constructed 
by  such  a  way  that  Eg  value  becomes  negligible  at  V  =  (3  +  5)  Vq 
[6]. 

The  temperatures  in  region  of  V  >  Vq  is 
T  =  T  +  E/C  . 

O  V 

Calculation  of  the  heat  capacity  as  a  function  of 
temperatures  Cv  =  f(T)  is  done  in  the  Debye's  approach 
[7]: 

Cv=  3  nk  [4  D(x)  -  3  x/(ex-1  )]  , 

x 

where  D(X)  =  — -  is  Debye's  integral,  X  =  9/T,  9  is 

X3  .  (e-1 ) 

o 

the  Debye's  temperature.  The  Debye's  integral  is  calculated  using 
standard  tables  with  linear  interpolations.  The  Debye's 
temperature  9  is  assumed  to  be  function  of  the  volume  only  and  is 
calculated  as  [ 8 ] : 
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|  ke  =  [  1  ,68/(Z  +  22)]  [p^Vd+p)2]  [ev], 

where  p  =  V  ,/V,  <p  =  0,6Z1/9,  V  ,  =  0,009  Z0,3/W, 

r  ref  1  ref 

Z  is  the  atom  number,  W  is  the  atom  weight. 

The  set  of  equations  (1)  -  (3)  is  solved  by  the  finite 

f  differences  method  on  the  Lagrangian  triangular  greed. 
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APPENDIX  B.  SIMPLE  CODES. 

A  series  of  simple  codes  for  personal  IBM-compatible  computers  has 
been  created  to  estimate  an  equation  of  state  of  composites  and  to 
test  constitutive  relations.  All  codes  have  been  written  in 
FORTRAN  and  presented  in  diskette. 

B.l.  MIXTURE  is  a  code  to  calculate  Hugoniots  of  mechanical 
mixtures  through  known  Hugoniots  of  components.  It  works  together 
with  the  file  "HUGS . DAT"  which  contains  Hugoniots  of  series  of 
materials.  You  can  include  any  additional  Hugoniots  into  HUGS . DAT 
file  using  the  same  format.  In  this  case  you  have  also  to  change 
the  number  in  first  position  of  the  "hugs.dat"  file  (amount  of 
Hugoniots  in  the  list) .  To  use  this  code  just  do  run  it  and  follow 
to  questions.  The  system  of  units  is  SI. 

B.2.  MIXSOUND  is  a  code  to  calculate  the  sound  velocity  as  a 
function  of  the  content  of  components  in  a  mixture.  It  works  with 
the  "HUGS.DAT"  file  too. 

B.3.  EOS  is  a  code  to  check  an  equation  of  state  which  was  used  in 
1-D  calculations.  It  is  the  equation  of  state  of  the  Mie-Gruneisen 
type.  The  cold  compression  isentrope  is  calculated  through  the 
Hugoniot  of  matter  basing  on  coinciding  of  the  Hugoniots  and 
isentrope  in  the  pressure  -  particle  velocity  coordinates: 

p  C2  7-7 

r0  0  0 

p  =  -  [exp(4S  - )  -1], 

4S  7 

o 

where  Cq  ,  3  are  coefficients  of  the  linear  relationship  between 
shock  velocity  and  particle  velocity.  The  Gruneisen  parameter  is 
assumed  to  be  constant.  This  simple  equation  of  state  gives  enough 
good  description  of  compressibility  of  solids  and  sound  velocity 
near  the  Hugoniot.  Results  of  calculations  with  EOS  show 
variations  of  the  Gruneisen  parameter  along  Hugoniot  when  this 
equation  of  state  is  used. 

B.4.  COMPOS  is  1-D  Lagrangian  code  used  to  check  constitutive 
relation  which  describes  a  wave  dispersion  in  composites.  It  works 
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together  with  COMPINST  file  where  all  parameters  of  calculations 
are  established.  COMPINST  contains  commentaries  to  all  parameters. 

B.5.  DUBRI  is  1-D  Lagrangian  code  to  calculate  shock-wave 
processes  in  the  system  impactor-target-window.  Three  different 
materials  can  be  included  in  calculations  and  each  of  them  can  be 
ductile  or  brittle.  This  code  works  together  with  the  file 
"PARAMS"  where  all  necessary  parameters  are  installed.  The  system 
of  units  is  SI.  Commentaries  to  parameters  are  given  in  the  PARAMS 
file.  The  elastic-plastic  properties  of  materials  are  described  by 
the  structural  Mar zing  model  which  has  been  presented  in  part  4.3 
of  this  report.  Two  parallel  elastic-plastic  elements  have  been 
incorporated  into  the  code. 
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Figure  1.  Sound  velocity  in  the  aluminum  and  tungsten  mixture 
calculated  in  additive  approach. 
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Figure  2.  Hugoniots  of  mechanical  mixtures  of  tungsten  with 
aluminum,  tungsten  with  epoxy,  and  aluminum  with  epoxy 
Volume  fraction  of  heavy  component  is  0.375. 
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Figure  3.  Modes  of  deformation  of  a  brittle  material. 
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Figure  4.  Deformation  of  the  composite  material  in  shock  wave 
Boundary  velocity  is  1  km/s. 
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Figure 


.  A  stress  wave  evolution  in  composite  at  100  ns  of 
the  initial  rise  time.  Hydrodynamic  approach. 
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Figure  6.  Comparison  of  the  shock-wave  stress  profiles  in  the 
composite  calculated  in  hydrodynamic  (points)  and 
elastic-plastic  (lines)  approaches. 
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Figure  7.  Shock  compression  of  porous  aluminum. 
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Figure  8.  Uniform  dynamic  compression  of  porous  tungsten  at  1 
km/s  of  boundary  velocity. 
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Figure  9.  Average  stresses  in  porous  tungsten  as  a  function  of 
the  strain  at  uniform  dynamic  compaction. 
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Figure  10.  Shock  compression  of  hard  ceramic  with  shear  bands. 
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Figure  11.  Average  stresses  in  the  hard  ceramic  as  a  function  of 
strain  at  the  uniform  dynamic  compression.  Friction 
coefficient  was  varied  from  0  to  0.33. 
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Figure  12.  Compression  wave  profiles  in  hard  ceramic  as  a  result 
of  2-D  simulation.  Solid  line  shows  the  elastic-plastic 
approach.  Dotted  lines  show  results  of  calculations 
with  intersecting  shear  bands.  In  case  of  the  dashed 
lines  the  shear  bands  were  not  intersected. 
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Figure  14.  Pressure  histories  in  the  middle  sections  of  two  layered 
composites. 
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Figure  15.  Results  of  the  pressure  profile  measurement  with 

manganin  gauge  at  the  interface  between  the  layered 


copper-polyethylene  composite  and  a  copper  plate. 
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Figure  16.  An  influence  of  stress  relaxation  in  polyethylene  of 
the  layered  copper-polyethylene  composite  on  the 
pressure  profile  at  the  interface  with  a  copper  barrier. 


Figure  17.  An  influence  of  elastic-plastic  properties  of  the  layered 

aluminum-tungsten  composite  on  the  pressure  profile  at 
interface  with  a  tungsten  barrier. 
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